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1.0 SUMMARY 


An improved panel method for the solution of three dimensional flow about wing 
and wing-body combinations with leading edge vortex separation is presented. 
The method employs a three-dimensional inviscid flow model in which the 
configuration, the rolled-up vortex sheets, and the wake are represented by 
quadratic doublet and linear source distributions. The strength of the 
singularity distribution as well as shape and position of the vortex spirals 
are computed in an iterative fashion starting with an assumed initial sheet 
geometry. The method calculates forces and moments as well as detail surface 
pressure distributions. Improvements include the implementation of improved 
panel numerics for the purpose of eliminating the highly non-linear effects of 
ring vortices around doublet panel edges, and the development of a least 
squares procedure for damping vortex sheet geometry update instabilities. 

The documentation is divided up into two parts: 

Volume I Theory Document 

Volume II User's Guide and Programmer's Document 


Volume I contains a complete description of the method. A variety of cases 
generated by the computer program implementing the method are presented. 

These cases are of two types. The first type consists of numerical studies, 
which verify the underlying mathematical assumptions of the method and 
moreover show that the results are strongly invariant with respect to such 
user dependent input as wing panel layout, initial sheet shape, sheet rollup, 
etc. The second type consists of cases run for the purpose of comparing 
computed results with experimental data, and these comparisons verify the 
underlying physical assumptions made by the method. 

Volume II contains instructions for the proper set up and input of a problem 
into the computer code. Program input formats and output are described. A 
description of the computer program and its overlay structure is also 
presented. 



2.0 INTRODUCTION 


2.1 Background 

The flow at the leading and tip edges of a swept wing with sharp edges 
separates at moderate to high angles of attack, the separation producing 
vortex sheets that roll up into strong vortices above the upper surface of 
the wing. The formation of these vortices is responsible for the 
well-known nonlinear aerodynamic characteristics exhibited over the 
angle-of-attack range, (Figure 1). Experimental studies (ref. 1) of the 
vortex sheet separating from a slender sharp-edged wing revealed that the 
rolled-up part of the vortex sheet consists of three regions: an outer, 

convection dominated region in which the distance between turns is large 
compared to the diffusion distance; and an inner region where the distance 
between turns is of the same order of magnitude as the diffusion distance; 
and inner, diffusion-dominated, viscous core which is very small, 
representing only about 5 percent of the vortex diameter. In addition, 
studies (refs. 1, 2) of the principal vortex indicate that its shape and 
strength are relatively independent of Reynolds number. The relative lack 
of viscosity dependence suggests that the flow may be regarded as 
potential, with the free shear layer represented either as a vortex sheet 
or, equivalently, a doublet distribution, supporting a discontinuity in 
tangential velocity. 

Many theoretical methods have exploited this fact to predict various flow 
characteristics. The leading-edge-suction analogy described in references 
3, 4, and 5 provides a method suitable for calculating the magnitude of 
the nonlinear vortex lift on a rather broad class of wing planforms. 
Polhamus (ref. 3) reasoned that the normal force needed for the flow 
around a leading edge to reattach to the wing is equivalent to the leading 
edge suction force necessary to force the flow to be attached to the 
leading edge in an unseparated condition. The unseparated leading edge 
suction force is calculated, and is then rotated normal to the wing to 
obtain the lift contribution of the leading edge vortex. The total wing 
lift computed by this method agrees well with experimental data, but the 
leading-edge-suction analogy does not give flow-field details or detailed 
surface pressure distributions. Several attempts had been made in the 
past toward the theoretical prediction of detailed pressure distributions 
and flow fields about swept wings with leading edge vortex separation. 

Most of these past methods are limited to slender configurations, a 
considerable simplification because the problem can be reduced to a 
solution of Laplace's equation in the cross flow plane, for which 
conformal mapping becomes a powerful tool. Smith (ref. 6) developed the 
best known method of this type by improving the work done earlier in 
collaboration with Mangier (ref. 7). Assuming conical flow, which is 
approximately valid near the apex of the wing, he was able to predict 
qualitatively the type of pressure distributions that had been observed 
experimentally. Those pressure distributions exhibit a vortex-induced 
pressure peak at about 70 percent of the local semi span of the wing. 

Toward the trailing edge. Smith’s method overpredicts the experimental 
load distribution by a considerable amount, because the conical theory 
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FIGURE 1 LEADING EDGE VORTEX FLOW 
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does not satisfy the Kutta condition at the trailing edge. Conical flow 
methods were followed by fully three-dimensional techniques in which the 
vortex is represented by single or multiple line vortices (refs. 8-13) or 
by a vortex sheet (refs. 14-16), and in which the trailing edge Kutta 
condition is enforced. These methods have enjoyed reasonable success in 
predicting overall configuration forces and moments and in some cases wing 
pressure distributions. A current review of various methods is presented 
in reference 17. 

The method presented in this report is basically that of reference 15. * 
This method was originally developed by the Boeing Company under contracts 
NASl-12185 and NASl-13833 from the Langley Research Center. The method is 
capable of predicting forces, moments and detailed surface pressures on 
wing and wing/body combinations assuming the separation lines are known. 
The wing geometry may be arbitrary in the sense that leading and trailing 
edges may be curved or kinked and the wing may have arbitrary camber and 
twist as long as in real flow it produces only a single well developed 
vortex system. The method employs an inviscid flow model in which the 
configuration surface is represented by source and/or doublet singularity 
panels, and the rolled-up vortex sheets and wakes are represented by 
doublet panels alone. The Kutta condition is imposed along all wing 
edges. Strengths of the singularities as well as shape and position of 
the free vortex sheet spirals are computed iteratively starting with an 
assumed initial sheet geometry. The original method had been in use for 
some time now with generally good results, however certain shortcomings 
had become apparent. First the iterations determining sheet shape and 
position became unstable in seemingly random cases, making parametric 
studies difficult (reference 18). Minor changes in wing paneling, for 
example, have sometimes caused a well converged case to diverge. 

Secondly, computed lift coefficients for wings of large aspect ratio 
tended to be higher than those predicted by the suction analogy, and 
experiment (reference 18), The effort to solve these problems became the 
focus of contracts NA$1-15169 and NASl-15275 from the Langley Research 
Center. A coordinated effort which also included work conducted for the 
Boeing Independent Research and Development Program was ultimately 
successful in overcoming these deficiencies. The effort to solve these 
problems is summarized in the following sections. For purposes of 
completeness the independent Boeing work is included in this documentation. 

2.2 Approach to the Problem 


The convergence problem was addressed first in the hope of creating a more 
reliable tool for investigating the aspect ratio problem. To improve 
confidence in the numerical features of the method a general upgrade of 
the numerics was made. The upgrade included such minor things as more 
precise calculation of the geometry and network edge matching 
sensitivities, but the major effort was the implementation of parametrized 
panels and doublet splines in order to ensure continuity of geometry and 
doublet singularity strength across all panel edges, thereby eliminating 
the highly non-linear effects of line vortices (discontinuities in doublet 
strength). This upgrade did indeed enlarge the class of problems over 
which convergence was achieved; nevertheless some rather simple cases 


4 



still diverged. It was therefore necessary to look at more fundamental 
possibilities. A detailed analysis of divergence indicated that because 
of certain paneling anomalies, satisfaction of some of these boundary 
conditions required rather substantial kinks in the vortex sheet locally 
which, as pointed out by Rubbert, set off a built-in instability in the 
vortex sheet updating procedure. A very simple least-squares penalty 
technique was developed to damp this instability with the result that 
convergence was achieved in all of a wide variety of previously diverged 
cases to which the technique was applied. 

It has been the author's belief that the lift coefficients calculated by 
the current method should tend to agree with those of the suction analogy 
wherever the assumptions of that theory are valid; hence the attack on the 
second problem began with studies designed to check the numerical 
implementation of the method. These studies included the determination of 
the effect of variations in panel density, panel layout, sheet roll-up, 
initial sheet shape, etc. In all cases the studies proved that the 
boundary value problem associated with the model was being solved quite 
accurately so that the model itself was in error. It was subsequently 
discovered that use of the linearized pressure formula in the wing wake 
(known to be somewhat inadequate at low aspect ratios) was causing 
substantial loss of the wing trailing edge Kutta condition at high aspect 
ratios. The use of a fixed design wake then eliminated the problem and 
produced excellent lift coefficient comparisons. 

The remainder of this report is organized as follows: Section 4 describes 

the general features of the current method as a point of departure for 
section 5 where the advances leading to the solution of the aforementioned 
problems are detailed. Section 6 gives examples of numerical verification 
of the method and Section 7 gives examples verifying the physical 
assumptions of the model. An appendix containing further details of some 
of the features of the method is also included. A users guide and 
programmers document is provided in a separate volume. 



3.0 NOMENCLATURE 


a 

A 

AR 

b 

B 

c 

Cl 

Cn 

Cm 

Cp 

D 

F 

F 

G 

H 

K 

A 

I 

L 

M 

oo 

A 

n 

->■ 

n 


tangent basis vector 
compressibility matrix 
aspect ratio 
local span 

boundary of fluid domain 
chord 

drag coefficient 
rolling moment coefficient 
lift coefficient 
normal force coefficient 
pitching moment coefficient 
pressure coefficient 
fluid domain 

equations determining singularity parameters 
force vector 

equations determining vortex geometry parameters 

hyperboloidal panel 

equations penalizing panel twist 

panel width 

unit vector along vortex core or network junction 

curve on B across which p is discontinuous 

number of grid point rows on a network 

free stream Mach number 

surface unit normal vector 

normal vector at panel center 
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NOMENCLATURE (CONTINUED) 


"c 

N 

P 

Pi 

P2 

-► 

P 

-► 

Q 

-► 

Qi 

Qo 

Qs,Qt ,0st 

R 

-► 

R 

s,t 

S 

-► 

V 

-► 

V 

-> 

Va 


Voo 

-► 

w 

-► 

w 

Wa 

A 

X 

x,y,z 

a 


surface co-normal center 

number of grid point columns on a network 

pressure 

isentropic pressure 

second-order pressure 

field point 

point on boundary B 

nine canonical panel points 

panel center 

parametric coefficients defining H 
compressible magnitude of R 
vector from Q to P 
hyperboloidal surface parameters 
singularity surface 
perturbation velocity 
total velocity 

average surface value of total velocity 
free stream velocity magnitude 
free stream velocity 
perturbation mass flux vector 
total mass flux vector 

average surface value of toal mass flux vector 
unit vector along x-axis 
Cartesian coordinates 
angle of attack 


NOMENCLATURE (CONTINUED) 


6 

Y 

Y ' 

A 

A 

6 

C 

n 

6 

© 

A 

A 

y 

(C , n,c ) 

p 

p 

Poo 

a 

Y 


/ 1 - M 2 

delta wing semi apex angle 
ratio of specific heats 

jump in quantity across singularity surface or line 
change in quantity from one iteration to the next 
flap angle 

surface vorticity vector 
span fraction 

vortex system orientation angles 

all vortex systems geometry parameters 

vortex system scale factor 

all singularity parameters 

doublet strength 

doublet strength at Q 

fed sheet scale factor 

chord fraction aft of trailing edge 

panel point coordinates in local panel 
coordinate system 

fluid density 

Newton iteration step size limiter 
free stream fluid density 
source strength 
perturbation potential 
gradient of perturbation potential 
gradient operator 
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NOMENCLATURE (CONCLUDED) 


co-gradient operator 
belongs to 

vector cross product 



4.0 DESCRIPTION OF THE METHOD 


4.1 Theoretical Model 


The essential elements of the present flow model, as outlined in Figure 2, 
are the configuration surfaces (wing, body, etc.), the trailing sheet 
(wake), the sheet emerging from the wing leading edge and tip (free 
sheet), and the rolled-up core or spiral region (fed sheet) fed by the 
leading-edge and tip-vortex sheets. The following boundary conditions are 
imposed on these elements: 

o The conf iguration surface must be impermeable. 

0 The free sheet and wake cannot support a pressure difference and must 

be impermeable as well. 

0 The fed sheet is an extension of the free sheet and feeds vorticity 
to the vortex core (modeled as a simple line vortex). The boundary 
condition governing fed sheet size and core orientation is that the 
total force induced on the fed sheet and core by the rest of the 
configuration be parallel to the core. The size of the fed sheet is 
chosen initially by experience or from the conical flow results of 
Smith (ref. 6). 

0 Kutta conditions are imposed along the appropriate leading, side, and 
trailing edges of the wing in the presence of free sheets emanating 
from these edges. 

4.2 Basic Concepts 


The Prandtl-Glauert equation 

+ + = ^ ‘^lo ( 1 ) 

is assumed to govern the perturbation velocity potential <p in the flow 
field about the configuration. Here ;yie x-axis is taken as the freestream 
direction, i.e. = Voo x, where Voo is the freestream velocity and v 

its magnitude. Total velocity V is then defined by v = Voo + where T 
= ( (0x> <^y’ *^z) ) is the perturbation velocity. The definition of 
impermeability and pressure appropriate to equation (1) is an open 
subject. The mathematically natural choice of zero normal mass flux and 
the second order pressure formula (reference 19) is preferred. The total 
mass flux vector W is defined as 

^=Poo^oo + W (2) 

where w is the perturbation mass flux vector defined by 

w = Poo i.0^ ^y’ ^ ^ ^ 
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• DIFFERENTIAL EQUATION 



FIGURE 2 FLOW MODEL 





To first order in perturbation quantities W is equal to pV (reference 
19) hence impermeabil ity can be expressed by 

(W-n) = 0 (4) 

where n is the surface normal. Equation (1), rewritten as 0, 

expresses conservation of mass, and equation (4) then guarantees that even 
if the configuration is such that the assumptions used to derive equation 
(1) are violated locally there is still no net production of fluid at the 
boundary surfaces. 


The second order pressure formula is 

~ Poo" [Poo(^oo ’"t) + V2(\ • w)J 

and it agrees with the isentropic formula 

M. 


Pi = Pc 


o V ^ 

r'oo ^ oo 

-yMoo^ 


(n-ti 


V 2 

^ oo 


2 7 

OO J' lx r1 n X r 2 
OO 


(IVi^-V 2 )]t-i -1 


(5) 


( 6 ) 


to first order in perturbation quantities. Mathematically the second 
order formula is closely associated with equations (1) an (4) in that 
equation (1) is simply the Euler-Lagrange equation for the Bateman 
variational principle, 

SSf P^ ~ stationary ^ ^ ^ 

for which specification of (iT • A) is the natural Neumann boundary 
condition. Of great importance in this case is that the second order 
pressure formula produces consistent force calculations for arbitrary 
configurations when force is defined in the usual way, i.e.. 


F=- Sf [V(W -(i) + p 2 n] dS 
S 


( 8 ) 


Equation (1) implies that F is zero when the surface S encloses fluid 
only, hence momentum is conserved exactly and the force on a given surface 
may be computed on any enclosing surface. 


Under rather general assumptions Green's third identify (references 20, 
21) shows that any solution of (1) at a field point P may be expressed as 
the potential induced by a combination of source singularities of 
strength a (Q) and doublet singularities of strength p (Q) on the 
boundary B of the fluid domain D: 

= + /; ^q(47r)<1So 

Here ^ is the position vector R is the compressible magnitude of R 
defined by 



A • ^ 


where A = 


( 10 ) 



0 0 \ 
0 ) 
0 /3V 


and ^ is the co-gradient with respect to Q defined by 

Vq 


( 11 ) 


The perturbation velocity v(P) associated with <f> may be computed by 
differentiating (9): 


t(P)= 4^a('5)Vj 

D 




dSQ + 


Xf 

B 



(12) 


whereupon application of Stokes' theorem to the second term on the right 
yields 

T(?)={^o(C>)Vp^5ijdSQ+ — 

^ 

Here ^(Q) is the surface vorticity vector defined by 


f = n © Vm 


(14) 


and L is any curve on B across which p has a discontinuity, say Aju. 

It is possible to show from equations (8), (9), (12), and (13) that across 
any singularity surface S 


A(p = fx 

(15) 

A 


n® A V = f 

(16) 

(n • AW) = Poo 

(17) 


(18) 


and 
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S L 


(19) 


where A denotes the difference between the value on the side on which n 
is defined and the other side, and the subscript A denotes the average of 
the two values. Here is the co-normal defined by 


= /52 A’^ n 


( 20 ) 


4.3 Mathematical Implementation of the Model 

The mathematical implementation of the boundary value problem set forth in 
section 4.1 is now described. Our basic unknowns are the source and 
doublet strengths on all surfaces and the position of the free, fed and 
wake sheet surfaces. To solve for these unknowns the following equations 
are derived. 


On a surface bounding the fluid on both sides (i.e., thin sheets such as 
the free sheet, wake and possibly wing) it is required that equation (4) 
hold on both sides so that from equation (17) it can be seen that such a 
surface is source free. Hence, these two boundary conditions can be 
replaced by the equivalent conditions that the surface be a doublet 
surface and that 

(WA-n) = 0 (21) 

On the wake and free sheet surfaces we have the additional requirement that 


A * ^''^A ® f)/(^‘n^) = 0 


( 22 ) 


Although both equations involve doublet strength and surface geometry, the 
primary function of equation (22) in conjunction with equation (12) is to 
define surface doublet strength whereas the function of equation (21) is 
to define the surface normal and hence surface geometry. An approximation 
often made by many methods for wake surfaces is that 

(23) 

in which case equations (22) and (21) determine wake vorticity and surface 
slope prior to solution. Such an approximation is not precisely valid but 
can nevertheless be made at least in the far wake because details of the 
wake flow there have little effect on wing pressures. However a more 
accurate representation is sometimes required in the near wake, primarily 
because the spanwise component of vorticity in the near wake (which can be 
large in that a portion of the wake is underneath the primary vortex core) 
has a strong influence on the Kutta condition at the wing trailing edge. 

In this case the solution approach used is to require that equation (22) 


14 



be satisfied on a fixed wake surfaces in the immediate vicinity of the 
wing trailing edge, which causes the wake vorticity to seek the correct 
lateral alignment there. 

The Kutta condition at the junction of a (thin) wing and vortex sheet can 
be stated in several ways; e.g., zero pressure jump at the wing edge, 
finite flow at the wing edge, no flow through the vortex sheet, etc. All 
of these phenomena are supposed to occur once the Kutta condition is 
satisfied, and which boundary condition is actually called the Kutta 
condition depends on which boundary conditions have previously been 
assigned to the wing and vortex sheet. In this case it has already been 
assumed that equation (21) holds on the wing and equation (22) on the 
sheet, neither of which guarantee finite flow at the wing edge. Infinite 
flow can be created only by a discontinuity in doublet strength or surface 
vorticity across the wing/sheet junction. A discontinuity in doublet 
strength creates a line vortex of strength equal to the discontinuity 
(equation 13). The powerful flow singularity induced by such a vortex is 
incompatible with wing impermeability and hence a discontinuity in doublet 
strength is already prevented by the wing boundary condition equation 
(21). However the weak fl^ singularity induced by a discontinuity in the 
surface vorticity vector f is not. In fact there exist vorticity 
distributions creating infinite velocities at the junction, yet no normal 
mass flux on the wing, nor pressure jump on the vortex sheet. These 
distributions are such that tends^to infinity as the junction is 

approached from the wing side, where i is the unit vector along the 
junction. Thus the choice as the Kutta condition (as originally proposed 
by Rubbert) is that f be continuous f^om the wing to the vortex sheet, 
i.e., that where A ^vortex - f wing . Other 

components of f may be discontinuous; however these discontinuities will 
be eliminated by the process of updating the vortex sheet to satisfy 
equation (21) (which can only happen if the surface normal n is continuous 
across the junction). 

Following Smith (reference 6), the purpose of the fed sheet is to condense 
the free sheet vorticity into a line vortex core, thereby terminating the 
free sheet rollup. Hence the fed sheet is chosen to be a doublet sheet 
whose strength is equal to the doublet value at the junction with the free 
sheet. Only the size and position of the fed sheet remain to be 
determined from boundary conditions. The boundary conditions are chosen 
to be consistent with those that would be applicable to the infinitely 
rolled up vortex sheet as well, namely^that the total force normal to the 
core be zero, i.e., 4®^ F = 0 where t is the unit vector along the line 
vortex core and AfiS given by equation (19). 

On surfaces bounded on one side by a non-fluid domain (e.g., body or thick 
wing) it is required that equation (4) hold on the fluid side. While this 
boundary condition formally completes the mathematical description of the 
boundary value problem set forth in section 4.1, it is not sufficient to 
guarantee a unique solution to the problem via equation (12) since both 
source and doublet strength on such surfaces cannot be determined by one 
boundary condition. However a (nearly) arbitrary boundary condition can 
be assigned to the other side (reference 22). In particular assigning 
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to 4 > the same value as on the fluid side implies p = 0 in view of 
equation (15) and leads to pure source surface modeling. Here, however, a 
doublet lifting surface is required inside the wing and a doublet lift 
carry-over surface is required inside the body. (Models using doublet 
singularities on the thick wing and body surfaces are also possible 
(reference 22) but have not yet been implemented). 

Summarizing the mathematical description of the boundary value problem 
described in section 4.1, the following equations determine singularity 
strength: 

A 

" *w = 0 on wing and body, p = 0 on body, a = 0 on thin wing 

ff = 0, =0 on free sheet and near wake 


Af-t=0 


on wing edges (Kutta condition) 


(24) 


a = 0, continuation of p on far wake, fed sheet, 
body wake, carry-over sheet 

Free and fed sheet geometry are then determined by 

A 

n *w^ = 0 free sheet 

£ ©AF = o fed sheet 

4.4 Numerical Procedure 


Solution of the boundary value problem of section 4.3 via equation (12) is 
accomplished with the basic panel method of reference 23. The method 
proceeds by dividing the boundary surface into networks. A network is 
defined as a smooth portion of the boundary which has subsequently been 
divided into panels and on which source and/or doublet splines have been 
defined accompanied by properly posed boundary conditions. The networks 
are assumed to be logically independent in that each network contributes 
as many equations as unknowns to the overall boundary value problem; hence 
networks can be added or dropped without total reformulation of the 
problem. Essential features of the computational scheme are summarized 
below. 

0 Geometry input for a network consists of a rectangular array of 

corner point coordinates. The portion of the surface lying between 
four adjacent corner points is approximated by an analytically 
defined panel. 

0 Discrete values of singularity strength are assigned to certain 

standard points on each network. These values are interpolated by 
source and doublet splines which on each panel are assumed to be 
defined by linear and quadratic distributions respecti vely. 
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o Certain standard points on each network are assigned as control 
points at which boundary conditions are applied. These points 
include panel center points as well as edge abutment points in the 
case of doublet networks. The latter serve to impose standard 
aerodynamic edge conditions automatically (for example, the Kutta 
condition, zero potential jump at thin edges, and continuity of 
singularity strength across abutting networks). 

0 The induced potential and velocity integrals of equation (12) 

(influence coefficients) are all evaluated in closed form, although 
standard far field expansions are employed when the control point is 
sufficiently distant from the influencing panel. 

Figure 3 displays the location of discrete singularity parameters and 
control points for various network types. These locations are selected to 
achieve singularity spline stability with respect to the type of boundary 
conditions applied at the control points (reference 24). Additional 
details may be found in appendix E. Figure 4 shows a typical 
thin-wing/body configuration paneling and Figure 5 shows the same 
configuration disassembled into networks. Control points located at the 
junction of two doublet networks (or at the junction of one network with 
empty space) are assigned to match singularity strength across the 
junction. If only one control point exists, doublet value is matched. If 
there are two opposing control points the component of vorticity along the 
junction is also matched. Control points at panel centers have the 
boundary conditions prescribed in section 4.3. 

In Figure 6 the free and fed sheet kinematics are illustrated. The fed 
sheet size and position in each x-cut are changed by varying the scale 
parameters A and v , the parameter A scaling the whole vortex system. 
The free sheet shape is changed by varying the panel orientation 
angles , keeping the relative lengths £j fixed. Note that the 
vortex system geometry has as many degrees of freedom as constraint 
equations (25). 


Let all geometry degrees of freedom be denoted by the vector & and all 
singularity parameters be denoted by A . The equation set (24) can now 
be denoted by 


F (A, ®) = 0 


(26) 


and the equation set (25) by 


G (A, ©1=^0 


(27) 


These equations are solved iteratively by Newton's method with controlled 
step size, i .e.. 
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FIGURE 3 LOCATION OF DISCRETE SINGULARITY PARAMETERS AND CONTROL POINTS FOR 
VARIOUS NETWORK TYPES 
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FIGURE 4 PANEL MODEL 
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( 28 ) 



where p represents symbolically the step size scaling parameter 6 
(see Appendix G). 6 is a positive number less than 1 and is chosen small 
enough to ensure a decrease in the norms of F and G. Note that the 
Jacobian matrix on the left requires differentiation of the panel 
influence coefficients with respect to changes in geometry. Details of 
this differentiation will be given in the appendices. 
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5.0 RECENT ADVANCES 


5.1 Improved Panel Numerics 

Equation (13) shows that a discontinuity in doublet strength produces a 
line vortex of strength equal to the discontinuity. The velocity flow 
field induced by such a vortex has a singularity proportional to the 
vortex strength and inversely proportional to the distance from the 
vortex. Moreover the sensitivity to a change in vortex position 
(appearing in the Jacobian of equation (28)) has a singularity inversely 
proportional to the square of the distance to the vortex. It is clear 
that if a control point is sufficiently close to such a vortex, the 
linearization implicit in equation (28) will be valid for small A® 
only, prolonging convergence. Fortunately the vortex core in the present 
model tends to stay a large (relative to its strength) distance away from 
wing and free sheet control points. However, our original panel 
discretization introduced unintentional line vortices which were of small 
strength, but which could be relatively close to wing and free sheet 
control points. These line vortices arose for two reasons. First, flat 
panels were generally used for efficiency reasons. These panels were the 
plane quadrilaterals formed by projecting the straight line segments 
joining the four corner points onto the plane passing through the 
midpoints of these line segments. This meant that gaps in geometry would 
be present at edges of panels belonging to networks where the corner 
points did not all lie in a plane, leading to ring vortices around each 
panel. Secondly our locally quadratic doublet distribution on each panel 
was defined by fitting a quadratic function to singularity parameters in 
an immediate neighborhood. Since each panel used different singularity 
parameters, continuity of doublet strength across panel edges could not be 
enforced, and this again led to line vortices. Even on fixed portions of 
the geometry these vortices created problems because the boundary 
condition at a panel center control point close to a discontinuity in 
doublet strength would occasionally attempt to suppress the singularity 
produced by this discontinuity, i.e., interact with the doublet spline to 
create greater continuity, rather than control finite flow in the 
appropriate manner. 

In order to alleviate this problem the hyperboloidal panels of Morino 
(reference 25) have been implemented along with a continuous doublet.^ ^ 
^line. The hyperboloidal panel H interpolating four corner points Q-, 
Q^, 0^ is defined as the point set 

H = (^(s,t) : 3(s,t) = se[-l,l],te[-l,l]} (29) 

where 

^ = ^'^(^l+Q2+Q3+d4) ^ = 

Qt = ^st = 


14($i-Q2-Q3+Q4) 

!4(5i-Q2'^33-^) 
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The panel H is depicted in Figure 7. 

^te_J;hat^H con^ins the straight line segments joining the edge midpoints 
Q^, Q -7 and the same as the flat panel, and hence H is also only a 
first order accurate approximation to the true surface^ H^ev^ H also 
contains the straight lines joining the corner points Q 4 

so it abuts adjacent panels exactly leaving no gaps. 

Doublet strength on H is defined parametrically in terms of nine doublet 
values Mi the points by the formula 

m(s, t) = M 1 [ */ 4 St( 1 +s)( 1+t)] + M2 l-s)( 1 +t) ] + M 3 [ !-«)( M)] { 30 } 

+Ma P/4St(l+s)(l-t)] +Ms [ViUl+tXlV)] +M6 H/2s(l;s)(l-t")] 

+M7 P/2t(l-t)(l-s^)l +Mg [‘/2s(l+s)(l-t")] +M5 [(l-s'')(l-t")] 


Note that along each of the line segments displayed in Figure 7 M(s,t) 
is quadratic and is determined solely by the (three) values of doublet 
strength at the midpoint and endpoints of the segment. The nine doublet 
values P-j of equation (30) are not all independent. On any given 
network the set of independent parameters determining the doublet spline 
on that network are the doublet values at the locations shown in figure 3, 
and in general these values do not include all the Pj of all the panels. 
Hence the nine doublet values y-j must be expressed in terms of the truly 
independent doublet singularity parameters. For this purpose an enriched 
set of network grid points (Figure 8 ) is defined which includes the 
original corner poiots, the panel edge midpoints and center points and 
therefore all the of al 1 the panels. At each of these points doublet 
strength is obtained by fitting a quadratic function to a sufficient 
number of neighboring singularity parameters by the method of weighted 
least squares. For stability the closest singularity parameters are 
weighted heavily, in particular doublet strength at an enriched grid point 
coinciding with a singularity parameter point is simply set equal to the 
value of that parameter. For enriched grid points along network edges the 
corresponding doublet values are allowed to depend only on singularity 
parameters located on that edge, and for this purpose a least squares fit 
based on arc-length along the edge is used. 

It is clear from the above construction that doublet strength will be 
continuous across panel edges in the interior of each network. At network 
junctions doublet strength can also be made continuous so long as the 
corner points and edge singularity parameter locations coincide. This is 
the case in figure 4 except for the wing/free sheet junction. 

(Modification of some of the edge singularity parameter locations for both 
doublet/design networks could eventually lead to precise continuity 
everywhere) . 
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Calculation of the influence coefficients for the above formulation is 
accomplished using equation (13), where the last term on the right may now 
be discarded entirely except for the fed sheet terminated edge. 

Evaluation of the remaining integrals is facilitated by an expansion which 
is similar to the curved panel expansion of reference 23, but which does 
not require the small curvature assumption. Rewriting the second integral 
on the right of equation (13) in terms of the parameters s and t we obtain 


H 


f(Q)®^Q(4lFR’‘*SQ 






® a to'- 




0 a 


so 


)1 dsdt 

(31) 


where 


^so = ^o-^t‘ 
^ to “ ^o~^s^ 


The term in square brackets on the right is simply a polynomial, but the 
integration cannot be performed in closed form since R is a quartic 
polynomial in s and t. However R can be approximated as follows. Let 
(s*,t*) minimize R on [-1,1] ® [-1,1], sojthat Q(s*,t*) is a closest point 
on H to P in the compressible norm. Let R*(s,t) be the quadratic part of 
R at the point (s*,t*), i.e., 

R*(s, t) = (Ro-QsS*-Qtt*~OstS*t*)-Qs(s-s*)-3i(t-t*^ (32) 


Then is continuous in s and t and may be approximated to any 

accuracy oy a polynomial T(s,t). Upon substitution of the approximation 


1 _ T(s,t) 

r 3 - R^r- 


(33) 


into the right side of equation (31) the integration may be performed in 
closed form. 
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Evaluation of the surface vortex integral using the above procedure is 
quite accurate but rather expensive. Hence a somewhat simpler 
approximation was developed for use in the intermediate field, i.e,, when 
the field point is a modest distance away from the panel but not 
sufficiently distant that a far field expansion converges. For this 
approximation the flat panel is used again along with the quadratic 
doublet distribution (based on surface coordinates ( 


M ,V) = Mo + ^ 




(34) 


Here the coefficients are obtained by expanding equation (30) in a 
Taylor's series about the panel center. (Note that the doublet strength 
obtained from equation (34) then agrees i identically with that of equation 
(30) along the lines s = 0 and t = 0). It was found that this 
approximation agreed well with the more exact calculation above - even for 
the field point at the panel center, hence it was decided to use it 
exclusively in the near field. While it would then seem that we are back 
to the flat panel/quadratic doublet distribution at least for the purpose 
of calculating influence coefficients it is important to realize that the 
underlying panel shape and doublet distribution are given by equations 
(29) and (30) respectively. With the old method it was impossible to 
ignore the line vortex term on the right side of equation (13) since the 
doublet strength could not be deduced from knowledge of the vorticity 
vector alone. Additional details on the hyperboloidal panels can be found 
in appendices A, B, C and D. 

5.2 Least Squares Geometry Update Procedures 

In Figure 9, a streamwise paneled delta wing with unwrapped free and fed 
sheets is shown. The control point on the first wing panel is displayed. 
The control points on the first row of panels on the free and fed sheets 
are also displayed along with their projection onto the wing. Since the 
first wing control point is so far from the apex, flow through the wing 
near the apex is not prevented and consequently the free and fed sheet 
control points in the first row encounter a somewhat different environment 
than those in subsequent rows. To satisfy the boundary condition at these 
control points the whole vortex system near the apex is required to move 
substantially inboard causing errors in wing pressures at the first wing 
control point. Obviously with such a paneling one cannot be too concerned 
with flow details at the apex anyway. In the past however, the flow 
anomaly there destroyed convergence everywhere even though the flow is 
better behaved farther aft on the wing. 

The problem was two-fold. First, all equations of the set (25) were 
required to be satisfied exactly. Because huge local anomalies in vortex 
sheet shape were required, the singularities themselves got heavily 
involved in solving equation (25) with consequential loss of stability. 
Secondly, the update procedure was such that local anomalies in sheet 
shape could propagate to other areas of the sheet; in other words the 
procedure itself was not fully stable. The basic reason has to do with 
the fact that for both the flat and hyperboloidal panels, the surface 
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normal at the panel center remains unchanged when the four panel corner 
points are alternati ngly perturbed equal distances above and below the 
average plane. For illustration purposes assume the free sheet lies in a 
plane as shown in Figure 10. The average mass flux vectors at the panel 
centers are depicted with arrows and are assumed to lie in the plane as 
well, except for the mass flux vector at the center of the first panel 
which because of a flow anomaly is assumed to be substantially out of 
plane. Assume that all corner points except for those in the inboard 
column (which are attached to the wing) may be perturbed in a direction 
normal ^o the surface (i.e., plane) and also assume that the mass flux 
vector W at each panel center^is little altered by such a perturbation. 

The boundary condition (w • n) = 0 is already satisfied on all panels 
except the first, where a perturbation of corner point 3 by a lar^ 
distance h is required to modify the panel center normal so that (W * ft) 

= 0. A perturbation of corner point 4 by -h is required to maintain the 
boundary condition on panel 2, and so on down the line. Thus the effect 
of the flow anomaly on panel 1 is propagated to the whole free sheet with 
the consequence that all the quadrilateral panels become considerably 
twisted. 


Probably the simplest method of damping with instability whenever it 
arises is to limit excessive panel twist. A measure of panel twist is the 
function 


K 


n*Q 


st 


(n-n) 


3/4 


(35) 


where n = Qs ® Ot* 
untwisted (flat) is 


The condition that all free sheet panels be 


K (® ) = 0 


(36) 


using the notation of equations (26) and (27). Equation (36) combined 
with equations (26) and (27) creates an overdetermined system of equations 
for A and 0 . View equation (26) as an equation which defines A as 
function of 0 , i .e. , 

F(A, 0 ) = 0 => A= f(0) (37) 


Substituting equation (37) into equation (27) results in having two 
competing equation sets for determining 0 , i.e.. 


G(f(0), 0) = O and K(0)=O (38) 

This system is solved in a least square sense after suitable normalization 
to account for dimensional differences as well as desired weighting. 
Obviously the penalty equation (36) should not be weighted too heavily 
since a free sheet made up entirely of untwisted panels cannot in general 
be a good approximation to a stream surface. Fortunately a small weight 
is all that is required. The instabilities produced by a local flow 
anomaly are severe enough that a very small penalty on panel twist forces 
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FIGURE 10 GEOMETRY UPDATE INSTABILITY 


relaxation of the boundary condition causing the local anomaly. 


The procedure for solving the overdetermined equation set is iterative as 
before. At the beginning of an iteration equation (26) is solved for A 
as a function of the current O using Newton's method with controlled 
step size, i .e.. 


4^. AA = -PF 
6A 


(39) 


where p represents symbolically a step size scaling parameter which is a 
positive number less than 1 and is chosen small enough to ensure a 
decrease in the norm of F (see equations G.9 and G.ll of Appendix G). 

Upon obtaining convergence a new estimate for © is calculated by solving 
the equation 

( dG AL 

6® + y© 

6k 
6 ® 

in a least square sense, where the Jacobian on the left is evaluated at 
the point A = f(©) as determined from (39) and 6J_ is 
calculated from 6© 

6f 6f 6F _ ^ 

6a 6® 60 

It is assumed here that G and K have been normalized appropriately. 

Again, p represents a step size scaling parameter which is a positive 
number less than 1 and is chosen small enough to ensure a decrease in the 
norms of G and K (see Appendix G). 



6.0 NUMERICAL VERIFICATION 


In this section cases generated by the computer program implementing the 
current method are presented. The purpose of these cases is to display those 
numerical characteristics of the method which are important for establishing 
confidence in the computed results. 

6 .1 Effect of Wing Panel Densit y 

In Figure 11, the effect of wing panel density on vortex system geometry, 
wing pressures, and forces and moments is shown. The primary effect is 
associated with spanwise density. In Figure 11a two conically paneled 
wings having 6 and 12 panels spanwise are illustrated. In both panelings 
the spacing is non-uniform with a concentration of panels where they are 
obviously most required, namely outboard under the vortex. The effect of 
an increase in density is to move the vortex system slightly outboard as 
shown in Figure 11b, and the consequence is somewhat higher pressures on 
the upper surface outboard of the vortex core as shown in Figure 11c. The 
lift coefficient is correspondingly higher (see also Figure 19b). Greater 
spanwise densities have been run but the incremental effect is negligible 
compared with that shown in Figure 11. 

6.2 Effect of Wing Panel Layout 

In Figure 12 computed results for the two panel layouts usually employed 
on delta wings, i.e., streamwise and conical are compared. Both layouts 
have 64 panels. The least square procedure described in section 5.2 was 
required to obtain convergence for the streamwise paneling. Wing 
pressures and force and moment coefficients are displayed in Figure 12b. 
Pressures for the streamwise paneling have been interpolated to the 
control point locations of the conical paneling for comparison purposes. 
The two cases were run about a year apart and in the intevening time a 
study on initial sheet shape was made, hence the free sheets have somewhat 
different panel spacing. However comparisons of pressure and force data 
still show excellent agreement. 

6.3 Effect of Vortex Sheet Rollup 


The cornerstone of the current method is, of course. Smith's device for 
terminating the free sheet rollup, namely the use of a fed sheet whose 
position and size are determined by the same overall force condition that 
would be applicable to an infinitely rolled up free sheet. The point at 
which the free sheet rollup should be terminated by a fed sheet depends 
upon the sensitivity of wing pressures to further rollup. This matter has 
been investigated in detail by Smith (ref. 6) under the assumption of 
conical flow, and the standard amount of rollup employed by the current 
method is based on his results. To verify the application to fully 
three-dimensional flow, delta wings of aspect ratio .25, 1.0 and 2.0 have 
been analyzed with an additional 1800 of rollup. Results at AR = 0.25 
are shown in Figure 13 and indicate a slight increase in lift (4 percent) 
with increased rollup. The effect of rollup is much less at the higher 
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FIGURE 11 EFFECT OF WING PANEL DENSITY 
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FIGURE 12 EFFECT OF WING PANEL LAYOUT 
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aspect ratios (1 percent at AR = 1.0 and less than 0.5 percent at 
AR = 2.0). Results for AR = 2.0 are shown in Figure 14. The 
conclusion is that the standard rollup is generally adequate for all 
models. The slight improvements available at low aspect ratios do not 
seem to warrant the added complexity and expense of additional rollup. 

6.4 Effect of Vortex System Kinematics 


The current standard vortex system kinematics shown in Figure 6 is, of 
course, only one of many possibilities. A good alternative is the 
kinematics of Smith (ref. 6) shown in Figure 15. Here, in contrast to the 
standard kinematics, angles are fixed and lengths are chosen as free 
parameters. Smith's kinematics has also been coded into the computer 
program implementing the current method with results typified in Figure 
16. Comparisons have been made for small and large free sheet rollups, 
for core locations inboard and in the vicinity of the leading edge, and 
for wings of small and large aspect ratio, all with similarly close 
results. If any difference has been noticed, it is that Smith's 
kinematics seems to converge somewhat faster than the current standard 
kinematics with corresponding savings in run cost. 

6.5. Kutta Condition 


In section 4.3 the manner in which the current method enforces the Kutta 
condition at a wing edge was described. It was pointed out that the 
equation formally assigned as the Kutta condition is actually a 
restriction on singularity strength rather than the flow, namely that the 
component of vorticity parallel to the wing edge be continuous onto the 
vortex sheet. This condition when combined with the standard flow 
boundary conditions assigned to the interior of the wing and free sheet 
then creates those properties at the edge commonly associated with the 
Kutta condition, i.e., finite flow, zero pressure jump, etc. It was also 
pointed out that finite flow requires continuity of all components of the 
vorticity vector across the edge, which can only happen when the wing and 
free sheet adjoin smoothly. This fact seems somewhat inconsistent with 
our converged vortex system geometry (e.g. Figure 11b) which usually 
displays a large discontinuity in surface shape at the wing-sheet 
junction. However the discontinuity in surface slope is qualitatively no 
different than any other on the vortex system and can be reduced (and in 
the limit eliminated) by dense paneling. In Figure 17 the effects of 
finer vortex sheet paneling at the wing junction are shown. Note that the 
discontinuity in surface slope at the junction is considerably reduced. 

To reduce the discontinuity in surface slope in the neighborhood of the 
junction to the point where the surface would appear smooth to plotting 
accuracy would require extremely fine paneling at enormous expense. 
Fortunately such paneling is not required unless precise details of the 
flow in the neighborhood of the junction are required for reasons other 
than establishing the Kutta condition. This is because the global effects 
of the Kutta condition are already accounted for by our particular 
implementation. Note from figure 16b that dense free sheet paneling near 
the junction has little effect on lift and moment coefficients, vortex 
core position and strength, and pressure distributions except at the 
junction. 
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b. WING PRESSURES 




45 








O 

< 

QC 




< 

a. 

CO 




Q 

LU 

Q 


O 

C_> 


UD 


LU 


cxl 

=D 

CJD 


47 









6.6 Effect of Initial Free and Fed Sheet Shape 

Aside from free sheet rollup and panel density the initial guess of sheet 
shape and size has little bearing on converged results (assuming, of 
course, the boundary value problem has a unique solution) although 
convergence itself will be affected. Figure 18 illustrates this point. 

In Figure 18b a converged pressure distribution and force and moment 
coefficients for an aspect ratio .25 delta wing at 20o angle of attack 
and no yaw are displayed. The initial sheet shape is considerably 
asymmetric, however the converged sheet shape is quite symmetric given the 
fact that rollup on left and right sides are slightly different. More 
importantly, the pressure distribution and force and moment coefficients 
are very nearly symmetric. Note the close agreement in values with the 
solution employing a plane of symmetry in Figure 13. 
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7.0 COMPARISONS WITH EXPERIMENT AND OTHER THEORIES 


In this section we present computed results for the purpose of examining the 
validity of the theoretical model described in section 2.1. Many such results 
have been reported previously in References 15, 18, 26. Results presented 
here will be those that are primarily new in nature, although some previously 
presented results will be repeated for the sake of completeness. Particular 
note should be taken of reference 26 in which the capabilities of the method 
in analyzing cambered wings are illustrated. 

7.1 Delta Wings 

7.1.1 L ift Coefficients as a Function of Aspect Ratio 

In section 2.2 the problem associated with the method overpredicting lift 
coefficient at high aspect ratio, and how this problem was eliminated by 
use of the more accurate near wake (type 6, Figure 3) at the wing 
trailing edge was discussed. Resolution of the problem is shown in Figure 
19. The dashed line in Figure 19b shows the lift coefficients calculated 
using a far wake (type 8) only. (Results using 60 wing panels would be 
somewhat higher at the high aspect ratios). Insertion of near wake (shown 
in Figure 19a) yields lift coefficients which are in substantially better 
agreement with the suction analogy and experiment (refs. 2, 27, 28, 29, 

30) at the higher aspect ratios. These lift coefficients are based on a 
near wake of 0.5 root wing chord in length divided into three rows of 
panels. (Note that the free/fed sheet vortex system must extend to the 
end of the near wake.) In order to study the sensitivity of the wing 
pressures to near wake length, two additional near wakes were tested. One 
was 1.5 root chords long and made up of 5 rows of panels (Fig. 19a). The 
other was 0.1 root chords and made up of 2 rows of panels. Wing pressures 
for all three wakes were practically identical indicating, as stated in 
section 4.3, that the primary function of this wake is to establish the 
proper wing trailing edge Kutta condition. 

The question of lack of agreement between the suction analogy, experiment 
and the current method at very low (0.5,) or high (.1.5) aspect ratios 
still remains. Examination of the results shows that for a very low 
(0.25) aspect ratio the experimental data lies nearly half way between the 
results of the current method (low) and suction analogy (high). However 
analysis of the experimental (reference 29) data reveals that at the 20o 
angle of attack shown, an asymetrical vortex has developed in the real 
flow. This is indicated by the rolling moment at zero yaw that developed 
at angles of attack greater than about 16o, jhe theoretical models 
assume symmetrical vortices, 0. M. Luckring (NASA LRC) has computed 
solutions at several angles of attack using the current method. His 
results presented in Figure 20 show excellent agreement between the method 
and experimental data at the lower angles of attack. Beginning with the 
150 angle of attack the results start to deviate from the experimental 
data. Since this is approximately where the leading edge vortices become 
asymmetrical while the theory vortices remain symmetric, it is conceivable 
that this phenomenon could somehow be responsible for the deviation. 
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FIGURE 20 AERODYNAMIC CHARACTERISTICS OF ASPECT 
RATIO 0.25 WING 


Suction analogy results are also shown in Figure 20. These results also 
deviate from experiment at the higher angles of attack but in the opposite 
manner. E. C. Polhamus (NASA LRC) has suggested that the suction analogy 
may be high for these cases because the analogy assumes complete 
reattachment of the leading edge vortex and therefore complete recovery of 
the suction force. As previously shown in Figure 13 for low aspect ratio 
deltas at high angles of attack, the size of the free sheet which 
represents the vortex shear surface has become so large that it may be 
preventing complete reattachment and therefore complete suction recovery.. 

At the higher aspect ratios, the current method is in good agreement with 
the suction analogy but both methods predict higher lift coefficients than 
shown by the data. Examination of the experimental results shows a loss 
in lift at the higher angles of attack due to vortex bursting. Since 
neither the current method nor the suction analogy can account for this 
phenonema exact correlation with experimental data is not possible. 

7.1.2 Pressure Distributions 

One of the most important features of the method is its ability to compute 
surface pressure distributions. This capability sets it apart from most 
other leading edge vortex methods. A comparison of detailed pressure 
distributions are shown in Figures 21 and 22. These results were 
calculated with the earlier version of the method (reference 14). Lifting 
pressures, shown in Figure 21, on an aspect ratio 1.0 wing are compared 
with the experimental data of Peckham (reference 2). Although only 25 
wing panels were used on one half of the configuration, the completely 
three-dimensional non-conical load distribution was predicted well 
including the location of the vortex induced pressure peak and the 
decrease of the load toward the trailing edge. Figure 21 illustrates the 
method compared to the experiment of Marsden (Reference 32) for an aspect 
ratio 1.4559 delta wing at an angle of attack of 14o, jhe general 
agreement between the predicted and the measured pressures is quite good. 
The experimental results clearly show the effect of the secondary vortex 
separation, which takes place on the upper surface just slightly outboard 
of the main vortex. The presence of the secondary vortex raises the 
suctions near the leading edge and lowers the suction peak due to the 
primary vortex. The theoretical method does not model secondary vortex 
separation and, consequently, produces a slightly different shape for the 
pressure peaks. This type of discrepancy will be found in most 
test-theory comparisons. 

The method is also capable of calculating pressures on a wing at yaw as 
well as at angle of attack. Such calculations require that the complete 
wing be paneled and that no symmetry conditions be imposed. Figure 23 
shows a comparison of calculated lifting pressures with the experimental 
data of Harvey (Reference 32 ) on a delta wing with 80 leading edge 
sweep. The theoretical results clearly predict the asymmetric character 
of the lifting pressure distribution. The descrepancy due to the presence 
of the secondary vortex is quite evident in this comparison. Additional 
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FIGURE 21 LOAD DISTRIBUTION OF DELTA WING, a 20.0° 
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FIGURE 23 YAWED DELTA WING 




comparisons were made for increasing yaw angles until one of the leading 
edges became parallel to the freestream. Agreement with test data 
comparable to that seen in Figure 23 was found in every case. 

In order to validate the accuracy of the new code in reflecting 
aerodynamic differences due to camber changes, two wings first analyzed by 
Kuhlman (Ref. 18) were reanalyzed. These were delta wings of identical 
planform tested by Wentz (Ref. 27 ). The wing had an aspect ratio of 1.15 
and 740 sweep. One wing was flat while the other had conical camber. 

These wings are illustrated in Figure 24. For the calculations, the wings 
were considered to have zero thickness. Calculated and measured pressure 
distributions are compared for both wings in Figure 25. While the 
differences in the measured pressure distributions due to the camber on 
the upper surface is small, the greater differences in the experimental 
data between the flat and cambered wing are also reflected in the 
theoretical results. 

The code must further accurately predict drag increments between different 
camber wings. Drag pol ars are presented for the flat and cambered wings 
in Figure 26. Both calculated and measured results are shown. This 
comparison shows the drag differences to be predicted reasonably well by 
the improved code. At C|_ = 1.2 the predicted drag reduction for the 
cambered wing is 5.7 percent compared to 7 percent given by Wentz for the 
measured data. 

7.1.3 Wake Vorticity Rollup 

An interesting property of the solution produced by the current method was 
discovered with the use of the near wake. In real flow over a delta wing 
with a leading edge vortex, the wake behind the wing will roll up into a 
vortex rotating counter to that of the leading edge vortex. The 
phenomenon has been known for some years and was clearly evident in the 
experimental measurements presented at a recent AGARD symposium on high 
angle of attack flows by Hummel (ref. 33). Examination of the doublet 
strength in the near wake and connecting free sheet clearly reflects this 
behavior. Doublet strength is plotted versus span fraction on the wake 
and unrolled free sheet for several cuts behind the delta wing in Figure 
27. Near the trailing edge the doublet strength gradually rises toward 
the wing tip indicating an outward spanwise vorticity flow in a direction 
counter to that of the vortex core denoted by the sudden drop of doublet 
strength. (Note that the vorticity is the gradient of the doublet 
strength). Moving away from the trailing edge the variation of doublet 
strength becomes flatter except for the jumps behind the tip and vortex 
core. This indicates the concentration of vorticity into two counter 
rotating vortices just as in the real flow. Presumably it would be 
possible to replace the whole wake by such counter rotating vortices at 
about 0.5 root chords behind the wing in order to determine the effect on 
downstream components of a more complex configuration. 
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FIGURE 2H ASPECT RATIO 1.15 DELTA WING 
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FIGURE 25 WING SURFACE PRESSURE' DISTRIBUTION 


DRAG POLAR COMPARISON 



FIGURE 26 DRAG POLAR COMPARISON 









DOUBLET STRENGTH IN THE NEAR WAKE 
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FIGURE 27 DOUBLET STRENGTH IN THE NEAR WAKE 



7.2 Rectangular Wings 


In Figure 28a converged geometry for an aspect ratio 0.5 rectangular wing 
with separated flow around the side edges is illustrated. Corresponding 
pressure distributions ( ACp) are plotted in Figure 28b along with 
similar distributions for the same wing in an unseparated condition (which 
were generated by simply removing the vortex networks). Force and moment 
coefficients for rectangular wings of varying aspect ratio are shown in 
Figure 28c and comparisons with the suction analogy and experiment 
(reference 35) are good. Convergence for these cases required use of the 
least squares technique described in section 5.2. 

Spanwise core locations at the trailing edge vary from 84 percent semi span 
at AR = 0.2 to 94 percent semi span at aR =1.0 compared with around 70 
percent semi span for delta wings. Core vertical displacements at the 
trailing edge vary from 13 percent chord at AR = 0.2 to 18 percent 
chord at AR = 1,0. This compares with 4 percent root chord at 
=0.2 to 9 percent root chord at AR =1.0 for delta wings. Trailing 
edge vortex core strengths for the rectangular wings are 12 percent higher 
at AR = 1.0 and 54 percent higher at AR = 0.2 than those of delta 
wings. The net result is that the spanwise flow at the trailing edge is 
markedly lower for rectangular wings, making use of near wake unnecessary 
for them. 

7.3 Arrow Wings 

Solutions for an arrow wing configuration have previously been calculated 
using the old code, reference 15. These results which are still 
considered valid are shown here for completeness. The experimental data 
is for an arrow wing-body configuration by Manro, references 35 and 36. 

An attempt with the old code to analyze the complete wing-body showed 
unacceptable slow convergence. Instead the configuration was modeled as a 
simple wing as shown in the inset on Figure 29. A near wake (type 6, 
Figure 3) was used to insure the proper Kutta condition at the trailing 
edge. The comparison with experimental data for the flat wing 
configuration shows generally good agreement considering the crudeness of 
the theoretical model. 

An analysis was also made of the configuration with a trailing edge flap 
deflection. The lifting pressure comparisons are shown in Figure 30 which 
also includes results from an attached flow solution. While some 
discrepancies do exist, the LEV code results are in substantially better 
agreement with the experimental data than are those of the attached flow 
theory. Part of the descrepancy is obviously due to the presence of a 
secondary vortex. Near the wing tip, the simplification of the wing 
planform to a pointed tip instead of the actual clipped tip may also 
account for the poor test-theory correlation in that region. 
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. WING PRESSURE DISTRIBUTIONS (ACp) AT AR = 0.5 

FIGURE 28 CONTINUED 
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FIGURE 29 LIFTING PRESSURE DISTRIBUTION ON FLAT ARROW WING 






8.0 CONCLUDING REMARKS 


The advances described in this report clearly improve the usefulness of the 
current method in the study of separated vortex flow. In addition, the 
numerical examples set forth in section 6 give reasonable assurance that the 
data computed by the method faithfully reflects the underlying flow model. 

For those cases in which the real flow deviates significantly from the single 
well -formed vortex assumed in the flow model, the method will generally fail 
to converge. Finally, the results accumulated to date show that the flow 
model itself is representative of the physics in a wide variety of cases. 

However, much work is still needed to improve the method to the point at which 
it exhibits reliability comparable to that of an attached potential method. 
Despite the advances made, difficulties may still occasionally be encountered 
with seemingly well posed problems. 


Boeing Aerospace Company 
P. 0. Box 3999 
Seattle, Washington 98124 
July 1979 
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Appendix A 

Hyperboloidal Panel Geometry 


Definition 

Let Qj, Q2, Q3 and Q4 be four panel corner points as shown in Figure 

A.l. The hyperboloidal panel H interpolating these points is defined as the 

point set 


H 




t) 


'(s,t) = ^o ^ V ^st"^ * " ^ 


[-1,1], t G 


1 ) 


where 

=1/4 (?i + ^2 ?4) 

Qs = 1/4 (^1 - ■q 2 - Q3 + ^ 4 ) 

Qt = 1/4 (Qi + ^2 - Q3 - Q 4 ) 

Qst = 1/4 (Qi - Q 2 + ^3 - ^ 4 ) 


General Characteristics 

The fact that H interpolaj^s the comer_^oints is easily ^monstrated by 
^ting that Q(l,l) = Q{-1,1) = Q2, Q(-l,-l) = ^3 and Q(l,-1) = 

Moreover the boundary of H simply consists of the straight lines 
joining yie corner points, a fact that can be checked by noting that for fixed 
s = Sq, Oiso,t) is linear in t and hence is a straight line segment. 

Similarly for fixed t = to- It is clear then that there will be no gaps 
between panels. 

0 is the average of t^ corner points and lies on H since Q( 0 , 0 ) = Qq. 

U the vector from to the ^dpoint of the line segment joining 

to ^ and is the vector from ^ to the midpmnt ^ the line segment 
■pining to The line segments Q( 0 ,t) = % ^ Qt t 

s, 0 ) = + QgS, s e[-l,l] belonging to H are then simply the line 

segments joining the midpoints of opposite sides. 

H is flat if and only if Q^t Ijgs in Xhe plane of Qg and Qt, i.e.^the 
near plane. If_jin pa^icular = ^5 thpn H i^a tri^gle with Q3 = 

Similarly Qst = ^ 0 ^ 5.^3* Qst = -Qs<=^ Ql = 

Q2 and Qst = -Qt < — !^ Ql = Q4* Q$t = ® H is a plane 
parallelogram, and vice versa. 


Normal 


Let us define 


3 s 


Q + 


Qjtt = [Q(l.t) - Q(-l,t)] (A. 2 ) 
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and 


and 


and 


9Q(s,t) 

at 




^ [ QU.I) - Q(s.-d] 


n = = Q^.© + q^® q^^s + Qjt® V 

= Qj® Qt + (QjS - V^®^st 

t _ ®s®*t 


(A.3) 


(A.4) 


A 

n = 


"M " l^s®^ 


(A.5) 


^ ^ls,t) .is the unit _^pper surface normal to H at (s,t). We note that 
s ® Qt / ^ ® Qt simply the unit normal to the 


Then n = ru 
ft(O.O) = Ts 
near plane. 


Area Element 


The area element dA on H is given by 

dA = 1^1 dsdt = 1 i^0*a^l dsdt 


(A.6) 


Note then that 


ndA = rTdsdt = ^^0 Jj. dsdt 


(A.7) 


Surface Derivatives 

Let f(s,t) be defined on H, and let 7 be the gradient operator in global 
coordinates. Then 


and 


n07f 


1 



9f 

3s 



3t 


(A.8) 


n07fdA = 



it 

3S 


r — 

‘s 9t 


j 


dsdt 


(A.9) 
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Variational Properties 


In this section we compute derivatives^ of_^varigus q^ntities _^ith respect to 
panel deformations, i.e., changes in Q]^, Q2, Q3 or Q4. Let Q4 be 
one of these points. Then 




9Q(s,t) 

=Ife.+£.s + 

1 01 '*S1 

Si* ^ 

Sti^* ) 

(A.IO) 

where I 

is the 3x3 

identi fy 

matrix and 
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1 
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^01 " 
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S2 

4 S3 “ 
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S4 
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1 

1 


1 

Si " 

4 

S2 = 

- T- S3 = 

“ 4 

^s4 

4 


1 


1 

1 


1 

Si " 

4 

S2 

4 S3 

” 4 

S4 

■ 4 


1 


1 

1 


1 

Stl 

4 

St2 

- 4 St3 ' 

4 

St4 

= ■ 4 


Next we have 


3a, 


9Q,. 


= I (e . + e 4.. t) 

'si sti ' 


(A.ll) 


and 


d a. 


3^, 


(A.12) 


and 


= Q, (eti " (Si " Sti^’ " St (Si^ " Si*) (A.13) 

3q ^ 

where Qc, and are matrices defined as follows. Let v be the 
vector tvi, V2, V3). Then 


0 -S s 


s ° -''1 


-s ''i ° 


(A.14) 


From other considerations we have 


8tni AT 

rc - n 


3Q. 


an 

a^- 


and so 


an 

a?. 


4 [i - 

ini J ^Q. 


Applications to equations (A. 6) and (A.7) are obvious, i.e. 

a .. aT 9n 


dA = 


30 


30 ,. 


dsdt 


\r (ftdA) . 

3Qi 30 , 


dsdt 


Finally we deal with equations (A. 8) and (A. 9) differentiating (A 8) 
respect to using (A.ll), (A.12) and (A.15) we get y v ; 




3f 


as 

3f 

at 


3at 
3 0 , 
3?. 


i"i 3?; 


-^2 

|nr 




-4^ 




a n 
a.n 


a?. 


Differentiating (A. 9) with respect to "q. we get 


a 


(n0VfdA) 



(e,., + e, s) 


’ti 


3f / . 


dsdt 


(A.15) 


(A.16) 


(A.17) 


(A. 18) 
with 


(A. 19) 


(A. 20) 
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Appendix B 


Panel Singularity Distributions 


Doublet Distribution 


We define 9 canonical panel locations as shown in figure B.l. We assume that 
doublet strength at these points is determined by a QxN^j matrix B^j which 
relates these strengths to a vector of neighboring doublet parameters. 

Let 




(B.l) 


where 


P = ( ^2’ ^3* 


Then the distribution of doublet strength |j(s,t) on the panel H is defined by 

U(s,t) = [“T St (l+s)(l+t)] + Uj St (l-s)(l+t)] 

+ Mj[.^st (l-s)(l-t)]' +y^ [^- St (l+s)(l-t)] 

+ Mg [-J- t (l+t)(l-s^)] + Mg[^ s(l-s)d-t^)] 


+ M 


7 2 

+ Mg [(l-S^)d-t^)] 


t (l-t)(l-s^)] + Mg' -^[s(l+s)(l-t^)] 


(B.2) 


From equation (B.2) we obtain 

9M(s.t) 

as 

1*1 [4“ (l+2s)t(l+t)] (l-2s)td+t)] + Mg -^[(l-2s)t(l-t)] 

+ '^4[- “T (l+2s)t(l-t)] + Mg [- st(ln)] + Mg [- 4“ d-2s)d-t2)l 

B . 3 ) 

+ M; [st(l-t)] + Mg[-|- d+2s)(l-t^)] + Mg [-2s(l-t^)] 
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and 


"^1 [“4“ (^■^2t)s(l+s^+M2 Cl+2t)s(l-s)] + (l- 2 t)s(l-i 



- -g- [(l-2t)s(l+s)] + Mj 

[-^ (l+2t)(l-s^)] +Pg [ts(l-s)] 

(B.4) 



j^-ts(l+s)] 

+ Pg [-2t(l-S^)] 


Let 

q (s,t) be the vector of 

coe ff icients 

of 15“ in (B.2), then 



M(s,t) 

= T (s,t) • 

M 

(B.5) 



= 9q(s,t) 

9s 

- * P 

(B.6) 


at 

- a^(s,t) 

at 

• y 

(B.7) 

Differentiation of Doublet Distribution 



Let X be the full vector of singularity parameters. Then 3^d/9X is 
simply a matrix containing almost all zeros except for those columns 
corresponding to the singularities of which have a 1 in the appropriate 

row. Then 


9u(s,t) 

dt 

= ■q(s,t)Bg 

at 

(B.8) 


3M(s,t) 
3s dt 

3T(s,t) „ 
as “d 

ax 



(B.9) 
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(B.IO) 


3P(s,t) , 3q(s,t) g 

at aX 3^^ 3x 


Assume ^he variation in 7 has been calculated with respect to any given corner 
point Qi . We have 

=t(s.t)'f 

9Q^- 3 (B.ll) 


9^P(s,t) _ 3?(s,t)^ 

3s aJj 


3^P(s t) ^ 3q(s.t)^ 3M 

at a^,- 


(B.13) 


Source Distribution 

We assume that the distribution of source strength on the panel a (s,t) times 
the area Jacobian |a^0a^| is linear, i.e. 

CT(s,t) = i n(0,0) l ^ ^ j (B.14) 

|n(s.tTl 0 s t 

Let Bj be the matrix which relates the coefficient vect^ a = (Cq , as* 
a^) to an N$ vector of neighboring source parameters so that 

0= (B.15) 


To obtain Bg we construct 
transformation matrix App 



a local near plane coordinate system at Qq with 
such that 



(B.16) 
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Let TTijbe the locations o^ neighboring source singularity parameters in this local 
coordinate system with ,f||^ ,C|^ ). Then the coefficients Oq , 

0’s and 0"^ are obtained by least .squaring the distribution a(4»n) = 

+ 0^^ + a^ri to the points if with heavy weight on the point’' (0,0,0) . 

In the same manner as for equation (B.ll) we have 


where 

and 



(B.17) 


Calculation of the Derivative of Doublet Strength with Respect to the Panel 
Grid Points 


From equation B.ll, we have 


9v(s.t ) 


^(s,t) 



(B.18) 


Where f (s,t) is the doublet strength at a control point expressed in terms of 
the parameters s,t of a. hyperbolic paraboloid_^element , q(s,t) is the 

coefficient vector of f , and = BjX with |j = (pi, p 2 > Pg) 

corresponding to values of doublet strength at the 9^canonical points of H-P 
panel . 


Finite difference approximation will be used to evaluate the derivatives 
9p Q-j • *^or each panel of a doublet/design network with its geometry 
perturbed', an outer spline, giving the dependence of the doublet strength at 
the 9 canonical points of the panel on the 16 neighboring points with the 
specified singularity parameters (see scheme to construct the outer spline for 
doublet/design type I), will be constructed. The 9 doublet strength 
coefficients (ff = (p^ , P2,-..» pg)) for the panel can be obtained from 
either multiplying the matrix B by the vector consistipg of the values of 
the 16 singularity parameters or the direct solution A p“= X (i.e. B^j is 
the generalized inverse of A). _^Then we perturb the 16 neighboring points one 
at a time and calculate a new yp for each perturbed point. Using ^ and the 
associated for each pertu^ed grid point, we compute the finite 
difference approximation of Bp/SQi . For a given panel, the derivatives 

of the doublet strength with respect to the network grid points can be 
assembled as a 9x16 matrix. 










JQ, 






JQl 



^ 

^ 




(B.19) 
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Appendix C 

Control Point Locations 


Panel Center Control Points 


We assume that all panel center control points are located at (s,t) = (0,0). 
Network Edge and Corner Control Points 


The hypothetical location of all panel edge and corner control points are also 
defined in terms of s and t. For a control point located at panel corner Qj 
for examp^ we have (s,t) = (1,1). For a control point located at edge 
midpoint ^ we have (s,t) = (0,1). For edge and corner control points with 
real boundary conditions the actual control point location will not coincide 
with the hypothetical location but will be withdrawn into the panel slightly. 
For example if the hypothetical location is (s,t) = (1,1) the actual location 
will be (s,t) = (1 - 5 , 1 - 6 ). Here 6 is say 0.1. If the hypothetical 

location is (s,t) = (0,1) then the actual location will be (s,t) = (0,1 -6 ), 

etc . 


83 



Appendix D 

Potential and Velocity Influence Coefficients 


Prel iminaries 

Let C be the compressibility direction unit vector and 6 ^ 
Define the matrix A by 

A = 6^1 + (1-a^) cc^ 


= 1 - m2 . 


(D.l) 


and the matrix B by 


B = I + (3^-1) 


(D.2) 


Let P be a field point and let Q be a point on H. Let 


R = P - Q 


(D.3) 


We now define R, the norm of R, by 

R AR 


(0.4) 


Next define 


f4) M „<1 

I2J M .»>1 


(0.5) 


The source potential (^>5 induced by a source distribution a on H is 
defined by 




dS. 


(D.6) 
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The doublet potential 4>d induced by a doublet distribution n on H is defined 

by 



(D.7) 


where 





= B 


d 


The source velocity induced by a source distribution G on H is defined 
by 

~ ^p^s " " TtT ^p “R“ ^ 

H 






(D.8) 


The doublet velocity Vj induced by a doublet distribution p on H is 
defined by 


V 


d 


1 

4tt 

47T 


= = JL 

'^p^d 4 tt 
ff 'n®7M)® Vn 

H ^ 

Xf(n® 7U)®- 

H 


H 

^ ^ 4^ J ’Q 4-® 


L 


.3 ^^0 


^2 p R 

4i^J 



(D.9) 



The line integral on the right of (D.9) cancels with the corresponding line 
integral for adjacent panels and may be ignored except on the panel edge 
corresponding to the vortex core. For the purpose of evaluating the remaining 
integrals on the right sides of (D.6), (D.7), (0.8), and (0.9) we approximate 
H by a flat quadrilateral panel Z. Z lies on ihe ^lane passing through ^ 
having a constant normal n in the direction Z is obtained by 

projecting the corner points of^H onto this plane in a manner to be described 
later. A local compressible (x,y',z') coordinate system may be defined on Z in 
the following way. 



where (x,y,z) are the coordinates of a point in the global coordinate system 
and (xQ,yp,ZQ) are the components of Qq. The transformation T is 
defined the column vectors comprising i.e., 



(D.ll) 


so that u, V, and w are the coordinate axes of the local coordinate system. 
The vector w is defined by 


"w*an (D.12) 


where 


1 

a a 


n =B 


A 

n 
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The vectors u and v lie in the plane of Z and are orthogonal. 


We have 




We now express the integrals of (D.6), (D.7), (D.8) and (D.9) in terms of 
local coordinates. We have 


dS - ~u(^~y d^dn * Oid^dn 


(D.14) 


= Vu ’ - x'f + in ' . y ')2 + ( . ^,,2 


(D.15) 


A ' 


n -R = z /3"^a 


(D.16) 


(n ®Vm) (x)Vr 




‘n’ 


0) 


1? - 


( 0.171 
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Hence 


^ = - — ff— '•e'dri' 

® NirJJ R 

E 


° ht\JJ 


(0.18) 


(D.19) 


Vs * T 


^ Q rr„C^-x',n'-y',c'-z') d^'dn' 

Ntt JJ ~3 


(0.20) 


■ ‘t- r)%yc « (..^ 1 ) 


The singularity distributions o and p can be expresses in terms of ^ and 
T] , however since a and y are defined in terms of s and t, and s and t are 
not polynomials in ^ and ri , the resultnt distributions would not be 
polynomials in C ^nd ri . Therefore we approximate a by a linear 
distribution and p by a quadratic distribution in ^ and ri . The 
approximation is defined by requiring that the derivatives of these 
approximations evaluated at^ = H= 0 agree with the exact derivatives of a 
and p to the respective order of the approximation. For this purpose we use 
the following formulas: 


H' 


^0 



-n 


II ' 

S 3t 


il 

an' 



II 

at 



^0 )|{ 

2^f I 

3t2 J 


9s s 


(D.22) 
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where 


ae'3W 


^0 3s ' ( Cst ' ?s ) 3t 


3^f 


^t5s) 


-€ 

3s3t 


n ^ 1 
S3s3t J 


3^1 

dr\2 


= A 


“ [Mst-^ Ve)fs^(^s5st-5s^c) 

+ f2 II 7 f f 4 . f2 1 

•' t „2 * e =«,-2 * 5 s ;;, J 


H 

dt 


3s' 

1 


(^s^t ■ CtOs ) 

- -^o[2Cte,n^, -(n,5, + 


iQ^.\,0) . TQ^ (et-Hj-O) = TQj ( ^3^, Q) = TQ^^ 

Now note from equation (30) that at s = t = 0, 


|y = -1/2 Pg* 1/2 p 8 


|H = 1/2 U - - 1/2 p 7 
at 5 ' 


3 y = 
3 s2 
3^ 


P. + P„ -21^ 


^ = 1/4 ■ 
as3t 


p- - p_ + p. _ p^ 


3s at 

2_L s u + p _ 2P 
3t2 ^ ' 


(D.23) 
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so that upon combining (D.22) and (0.23), 

may be expressed in terms of the nine canonical panel doublet values , 

and the result substituted into (D.19) and (D.21). For the source , , 

distribution of equations (D.18) and (D.20) we assumeCJ(^, ri)=Q + 
from the beginning and determine by the process of least 

squares (ref. 14). The resultant integrals are evaluated in the same manner 
as for (ref. 14). 

The line vortex term of equation (D.9) may be evaluated directly. Without 
loss of generality we assume the panel edge corresponds to t = -1 in figure 
A.l, Then 

R = - sR, 


where 


R- = P - Q7 • R * 1/2 { - Q, ) 


(D.24) 


also 


dH = Rg ds 


Hence 


-1 


p(s) ds 


(0.25) 


To evaluate f we note that 


M(s) - + M 1/2 (M^- 113)5 + 1/2 (M4+ M3 - 2M7 )s‘ 


(0.26) 


and 


R(s) =Ja + 2bs + cs' 


where 


-vr ->T -^T ->■ 

a = Rq A R^. b = -R^ A R^. c = r‘ A R^ 


Thus we need to compute H(l,3), H(2,3), H(3,3), where 

1 


H (M 


.K)-f 


M-1 


ds 


-1 


(0.27) 
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We obtain from integral tables 


where 


Then 


and 


where 



H (2,3) = -^[bH (1.3) + E (1.1)] 

H (3.3) = i [h( 1,D - aH(l,3) - 2bH(2,3)] 


H (1,1) = tL log 


\ 

Rl + Aj J 


l ^>0 




< 0 


and 


_Llog A'^l ~ ^ 


^2 > 0 , ■ 


I 


(cs + b) 


(D.28) 


(D.29) 


(D.30) 


0 (D.31) 
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Appendix E 


Schemes to Construct the Outer Spline for Doublet 
Analysis (2), Doublet Design Type I (4) and 
Type II (6), Doublet Wake Type I (8), Type II (14), 

Type III (16), Type IV (10) 

The outer spline gives the relation between the doublet strength at 9 
canonical points on a Hyperbolic Paraboloid element (panel) and the 
neighboring singularity parameters. For each canonical point, a local tangent 
plane coordinate system (^,n) is set up with the given canonical point as 
the origin. Then a quadratic function 

f(5.n) = ^ -r a3 5^ + a45n+ - 5 “ (E.l) 

is fitted (least squares) through the projections of some selected neighboring 
points where singularity parameters are specified. The coefficients aQ 
gives the desired relationship of the doublet strength at the given canonical 
point in terms of the neighboring singularity parameters. For the 9 canonical 
points on a H-P element, we have 9 such coefficients aQ's which define the 
outer spline for a given panel. 

In the following, we will discuss schemes to construct the outer spline for 
various types of doublet singularity networks. 

For a given panel of a network, the 9 canonical points consists of panel 
corner points, midpoint of panel edges, center of panel. We will show how the 
doublet strength at each of these locations depends on the neighboring 
singularity parameters 


NW 


NE 



Figure E.l The 9 Canonical Panel Points 
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(E2-1) For network edge grid points, the associated doublet strength is 
assumed to depend only on the network edge singularity parameters. 

The network edges are parametrized by arc length t. A least squares 
fit of the function f(t) = a + bt + ct2 through the four 
neighboring singularity parameters locations with two on each side of 
a given edge grid point yields the dependence of the doublet strength 
at the given point on the neighboring singularity parameters. 
Constraints are to be set on the two close singularity parameter 
locations. 

(E2-2) For network edge midpoints, the value of the specified singularity 
parameter defines the doublet strength at each midpoint. 

(E2-3) For each interior grid point, the quadratic function f ,n) is 
fitted (least squares) through the 12 neighboring points (its 
projections on the tangent plane with the grid point as origin) where 
singularity parameters are specified. A constrained least squares 
fit will be performed with constraints set at the 4 center points of 
these panels having the given grid point as their common corner 
point. It means that 4 equations corresponding to those 4 center 
points will be satisifed exactly in solving the least squares problem. 

e.g., Qp? shown in Fig. E.2, a quadratic function will be 
fitted (least squares) through the 12 neighboring points 
where ^31’ ^12 » ^22 » ^32 > ^42* 

^13 > ^^23 ^33» -^43 » ^24 » ^34 » 

specified with constraints at the 4 center points having A22> 

^32» ^23» ^33> associated singularity 

parameters. 

(E2-4) For each interior edge midpoint, the quadratic function f(C ,n ) is 
fitted through the 12 immediate neighboring points where the 
singularity parameters are specified. Again, a constrained least 
squares fit will be performed with constraints set at the center 
points of the two panels sharing an edge on which the given midpoints 
is located. 

e.g., N shown in Fig. E.2, a least squares quadratic function 
will be fitted thru 12 neighboring points where ^32* 

^42. ^52. .^33> ^43» ^53> ^34» ^44» 

'^54* '^35> "^45» 55 specified with 

constraints at the 2 center points with the associated 
singularity parameters X43, X44 also for W shown in Fig. 

E.2, the 12 neighboring points are locations where X23» 

\33y \43» ^33 » ^24 » ^34 » ^44 > ^54 » 

^25» X35, A45, A55, are specified. The points 

where X34 and X44 are specified will be the locations for 
setting the constraints. 
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(E2-5) For panel center points, the value of the specified singularity 
parameters defines the double strength at each center points. 

The following table gives the dependence of each canonical point of a 
given panel on the 21 neighboring singularity parameter (see the 
panel marked with C, N, S, W, E in Fig. E.2). 
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(E4-1) 

(E4-2) 


(E4-3) 

(E4-4) 


(E4-5) 


For network edge grid points, the value of the specified singularity 
parameters defines the doublet strength at each grid point. 

For network edge midpoints, the associated doublet strength is 
assumed to depend only on the network edge singularity parameters. 

The network edges are parametrized by the length t. A least squares 
fit of the function f(t) = a + bt + ct 2 through the four 
neighboring singularity parameters locations with two on each side of 
a given midpoint yields the dependence of the doublet strength at the 
given point in the neighboring singularity parameters. Constraints 
are to be set on the two close singularity parameters. 

e.g., for N^i in Fig. E.3, the four grid points where 
^31> ^41* specified will be used 

in the least squares fit with constraints set at grid points 
where X 21 X 3 ]^ are specified. 

For special case such as N^ shown in Fig. E.3, the corner point 
where Xix specified will be treated as two separate but 
identical points. Likewise for the midpoint , the corner point 
where X 5 ^ is specified will be treated as two separate but 
identical points. 

For interior grid points, the value of the specified singularity 
parameters defines the doublet strength at each grid point. 

For each interior edge midpoint, the quadratic function f(^ ,ri ) is 
fitted through the 12 immediate neighboring grid points where the 
singularity parameters are specified. A constrained least squares 
fit will be performed with constraints set at the two grid points 
which are vertices of the edge containing the given midpoint. 

e.g., for N in Fig. E.3, the grid points ^'“22* ^32» 

^ 42 » ^^ 52 » ^ 23 » ^ 33 * ^ 43 » ^ 53 ' ^ 24 » 

“^ 34 > ^ 44 » ^^54 specified will be used in the 

least squares fit with constraints set at grid points where 
^33 and X 43 are specified. 

Again for special cases such as W 21 shown in Fig. E.3, each of 
those three grid points where An, A 21 , A 3 X are specified 
will be treated as two separate but identical points. 

For each panel center point, the quadratic function ( ^ , n ) is fitted 
through the 16 neighboring grid points where the singularity 
parameters are specified. A constrained least squares fit will be 
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performed with constraints set at the four corner points of the panel 
of which the given center point belongs. 


e.g., for C 
Xp2» ^ 

*45 » 

fit with 

^33 » 


32 » 


in Fig. E.3, the 16 
A' 


42 


52 


neighboring grid 

^33 » 

^ 25 * 
used in 


points where 


^43 


^24» ^34» ^^44» A54, ^25, ''35, 

A 55 are specified will be used in the least squares 

constraints setting at four corner points where 


43 


‘34 


A 44 are specified. 


For special case such as shown in Fig. E.3, each of those grid 
points where A 21 , X 31 , X 22 > ^13 specified will be 

treated as two separate But identical points and the corner point 
where X]^]^ is specified will be treated as two separate but 
identical points. Similarly for C 21 shown in Fig. E.3, each of 
those grid points where Aji, X 2 l> X 3 jL, X 4 ;^ are 
specified will be treated as two separate but identical points. 


The following table gives the dependence of each canonical point of a 
given panel on the 16 neighboring singularity parameters (see the panel 
marked with C, N, S, W, E in Fig. E.3). 
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(E 6 - 1 ) For network edge grid points lying on the edge with singularity 

parameters specified at grid points, the doublet strength are simply 
the values of singularity parameters as for network edge grid points 
lying on the edge with singularity parameters specified at edge 
midpoints and two extreme corner points, the doublet strength are to 
be found by using the Doublet/Analysis approach (see (E 2 - 1 ). 

(E 6 - 2 ) For network edge midpoints lying on the edge with singularity 

parameters specified at grid points, the doublet strengths are to be 
found by using the Doublet/Design Type I approach (see (E 4 - 2 ). 


For network edge midpoints lying on the edge with singularity 
parameters specified at edge midpoints and two extreme corner points, 
the doublet strengths ar simply the value of singularity parameters. 

(E 6 ~ 3 ) For each interior grid point, the quadratic function (^,n) is 

fitted through the 12 neighboring points where singularity parameters 
are specified. A constrained least squares fit will be performed 
with constraints setting at the two edge midpoints which are on the 
two edges having the given grid point as their common vortex. 


e.q., for Q33 shpwn in Fig. E. 4 , the locations where 

^ 22 » ^ 32 » ^ 42 » ^ 23 » ^ 33 » ^ 43 > .^^ 24 ’ 

^34 » ^44 > A 25 » ^35 > “^45 specified will 

be used in the least squares fit with constraints set at edge 
midpoints where ^33» ^34 specified. 


For these interior grid points lying on the columns next to the first 
or last column, those network edge grid points will be used in the 
least squares fit. 


e.g., for Q22 shown in Fig. E. 4 , the edge grid points 
where ^11» ^31 specified will be used in 

conjunction with those edge midpoints where Xi2» 

^22 » ^32 » ./l3» ^23* ^33 » ^14 » ^24 » 

A34 are specified. 

(£ 6 - 4 ) For the interior midpoints lying on rows of the network, the doublet 
strengths are simply the values of the specified singularity 
parameters. 


For an interior midpoint lying on columns of the network the 
quadratic function f(^, n) is fitted through the 16 neighboring 
locations where the singularity parameters are specified. A 
constrained least squares fit will be performed with constraints 
setting at the four edge midpoints belonging to the two panels both 
having the given point as their common edge midpoints. 
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e.g., for N shown in Figure E.4, the locations where 

^ 22 » ^ 32 » ^42 » ^ 52 * ^ 23 * ^ 33 * ^ 43 » 

^ 53 » ^ 24 » ^ 34 » '^ 44 * ^ 54 » ^ 25 ^ ^ 35 » 

^ 45 * ^55 specified will be used in the least squares 

fit with constraints set at the four edge midpoints where 
^ 33 * ^ 43 » ^ 34 » ^44 specified. 

For those edge midpoints lying on the columns next to the first or 
last columns, the network edge grid points will be used in the least 
squares fit, e.g., for Np 2 shown in Figure E.4, the four edge grid 
points, where X^, ^? 1 » ^31» ^41 specified will 

be used in conjunction with those edge midpoints where ^ 12 » 

X 22 , ^32, ^42, X 13 , ^23* ^33* ^43» 

^14» ^24* ^34» ^^44 specified. In this case, the 

constraints will be set at the four edge midpoints where X 22 , 

^32 ♦ ^23 » ^33 specified. 

Special case such as the interior midpoints lying on the corner 
panels, each of the locations along the network edge where the 
singularity parameters are specified will be treated as two separate 
but identical points. 

e.g., for N ^^2 shown in Figure E.4, the locations where 

^ 12 » ^13» ^14 specified will be 

counted as twice in the formulation of the least squares problem. 

(E6-5) For each panel center point, the quadratic function f(^,n) is 
fitted through the 12 immediate neighboring locations where the 
singularity parameters are specified. A constrained least squares 
fit will be performed with constraints set on the two edge midpoints 

on the same panel that the given center point belongs to. 

e.g., for C in Figure E.4, the locations where ^23* 

^ 33 » ^ 43 » ^ 53 » ^ 24 » ^ 34 » ^ 44 » ^ 54 » 

'^25* '^55 specified will be used in 

the least squares fit with constraints setting at the two edge 
midponts where A 34 , X 44 are specified. 

For these center points lying on the first or last column panels, the 
network edge grid points will be used in the least squares fit. 

e.g., C 21 in Figure 4, the network edge grid points, where 
^ 11 » ^ 21 » ^31 > ^41 specified will be used in 

conjunction with the edge midpoints, where Xi 2 » ^ 22 » 

^ 32 » ^42 > ^ 3 » ^ 23 » ^ 33 » ^43 

specified. 

Again special case such as those center points lying on the corner 
panels, each of the locations along the network edge where the 
singularity parameters are specified will be treated as two separate 
but identical points. 
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e.g., for Cn in Figure E.4, the locations where An, 

^12* ^13» ^14 specified will be counted as twice 

in the formulation of the least squares problem. 

The following table gives the dependence of each canonical point of a 
given panel on the 20 neighboring singularity parameters (see the 
panel marked with C,N,W,S,E in Figure E.4). 
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Doublet/Wake Type I (8) 

The singularity parameters are specified at the locations as illustrated below. 



M 



Figure E.5 Doublet /Wake Type I (8) Network 
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First one computes the values of doublet strength at network edge grid points 
using Doublet/Analysis approach (see (E2-1). With the assumption of constant 
singularity strength along the streamwise direction, we define 

Mnw = uu = usw = ^C,= M.S “ 

and = ^SE = KQi+l) where M(Qii) and 

M(Q-}+1 i) are values of doublet strength at edge grid points Q-jj and 
Qi+1 i-’ The exceptions to the above definition are = )^W 

= = A for the first panel and nw ~ Mg = = 

7 for the last panel as illustrated in Figure E.5. 


Doublet/Wake Type II (14) 

The singularity parameters are specified at the locations as illustrated below, 



Figure E.6 Doublet/Wake Type II (14) Network 
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This type of network is used for fed sheet. Similar to Doublet/Design 
network, the values of the specified singularity parameters define the doublet 
strength at edge grid points. The values of doublet strength at edge mid 
points are computed using Doublet/Oesign approach (see (E4-2). With the 
assumption of constant singularity strength along the chordwise direction, 

(row M), we define = Xi. Mw = PC = 

^E = M(Mi) and = MS = ^^SE = X-j+i where A-j, Ai+i are 
values of singularity parameters specified at grid point. Q-j, Q-j+i and 
are edge midpoint on the edge joining and Q-j+i. 

Again using the assumption of the constant doublet strength along chordwise 
direction, the singularity distributions on the panel at 2nd row will be same 
as those on the panel at 1st row. 


Doublet/Wake Type III (16) 


The singularity parameters are specified at the locations as illustrated below 



Figure E.7 Doublet/Wake III (16) Network 
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This type of network is identical to Doublet/Wake Type II (14) and is used 
mainly for wake network attached to Doublet /Design Network. 


Same as Type II (14) network, the values of the specified singularity 
parameters define the doublet strength at edge grid points. The values of 
doublet strength at edge midpoints are computed using Doublet/Design approach 
(see (E4-2). With the assumption of constant singularity strength along the 
streamwise direction (row M), we define MnW = = MSW = 

A^, Mm = Uq = ~ and M|^£ = = 

^SE = Ai+i where A -j , ^i+i 3re values of singularity 
parameters specified at grid point. Q-j , ^*^9® midpoint 

on the edge joining Q-j and Qi+1. 


Doublet/Wake Type IV (10) 


There is only one singularity parameter specified at the location shown below. 



Figure E.8 Doublet/Wake Type (IV) (10) Network 
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The doublet strength is assumed to be constant over the whole network. The 
values of doublet strength at the 9 canonical points are all equal to the 
value of the specified singularity parameter A . 
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Appendix F 

Function Evaluation and Jacobian Formulation 
Function Evaluation 

The boundary conditions can be expressed symbolically in terms of the 
following equations: 


F(A»0,X,v) 


r 



0 

0 


- 0 


wing and body 

Free sheet and wake (F.l) 

Kutta condition 


G(A,0,A,v) ■ W - 0 
n 

H(A,0,X,v) ■ f * 0 


Free sheet 

(F.2) 

Fed sheet 

(F.3) 


Where A denotes the singularity (doublet and source) strength parameters, 0 
represents the angle of inclination of panel edges in transverse cuts defining 
the spatial location of free sheet, and A , v are free geometry 
parameters controlling the shape of free and fed sheets. 


The function F symbolizes the impermeability condition of wing and body, zero 
pressure jump of free sheet and wake, Kutta condition. The function G is the 
impermeability condition of free sheet. Finally the function H represents the 
global boundary condition of zero net force acting on the fed sheet and the 
line vortex. Each of these functions is evaluated in the following manner. 


1 . 


Wn -PooCV^ + BV^)-S 


at a given control point (doublet network) 


where 7oo is the free stream velocity Ya is the average perturbation 
velocity at the control point induced by al f panels, and n is the unit 
normal at the control point on the HP surface. The matrix B is defined by 


B =■ I + (0^ - 1) c c 0V^ (F-4) 

where 02 _ i is the compressibility direction unit 

vector and A = 3^1 + (i _ 3 2)^ ^T, por a control point that 
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lies on a source network to satisfy the upper surface impermeability, we 
have instead the following expression 


w„ - p<„(V„+ BV^)-n + JO P, 


(F.5) 


2. 


where a is the source singularity strength. 

ACp • -2(V„ + BV^)-(ft®7u) ® ri 


(F.6) 


where and B are defined as before, and n is evaluated as 

follows: 




(F.7) 


The unit upper surface normal on the HP surface is given by 

A as®*t 


1 ® a^I 


3. The force across the fed sheet is calculated from equation (19) with 
0=0. We require that the two components normal to the core be zero. 
In each column of fed sheet panels to be specified we require that 

+ Ju T^®dlj= 0 

A (F.8) 

\ ® tds + j" M w ^ X dtj = 0 

t _ L 


where S is the panel column surface, L is the line segment comprising the edge 
of the panel column (corresponding to the core) and the subscript e denotes 
quantises eva'^uated at the midpoint of L. Here is given by J ® 

wherej^ = £, jtis the vector line segment L. From*^equation (18) w# have^ 


“ X C) “ 


* ^ 
Apn 


( w .-n)(C X n) 


/ A 

(n 


n) 


(F.9) 


For modest Mach numbers n ®n ~ 0.; moreover by construction of the fed 
sheet doublet spline ~r is^ approximately perpendicular to the core so that 

A ^ 

Vg ® C a 0 • 
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Hence 


A 

- V. 




A A 
Vg-n 


and 


-ng (Wft® C )~*P "e" 


(F.IO) 


Thus (F.8) is approximated by 
fv = XO. 


1 A 

^ e 

e 




= 0 


= 0 


L Sheet 

E Panel Column | (F.ll) 

A ^ ' vr .ix; 

APi A, n, +U^w^®te 

Fed Sheet 
" Panel Column 

Here pj is evaluated at center of the ith panel of the panel column and A 
is the area of the ith panel. 


Jacobian Calculation 


Taking the partial derivatives of the functions F, 6 and H with respect to the 
variables A » 0 , A and v we have the Jacobian matrix 


V 




3W 

n 

\ 

3A 

30 

3X 

3v 

9C 

P 


as 

acp 

3A 

30 

3X 

3v 


5«nk 

3W . 
nk 

3W . 
nk 

3A 

30 

3X 

3v 

3f 

3f 

3f 

3f i 

, 3A 

30 

3X 

3v / 


(F.12) 


This gives the variations of the boundary conditions due to the perturbation 
of singularity strength and geometry parameters 

(1) 


3A 3A 



. n 


(F.13) 


no 



For body-source network, 3W^ 3[(V^ + + l/2o] 


3 A 2 3A 


(F.14) 


. 0 . 6,x,v 

30 .30 

3V, ^ aft 

. B —4. .ft + (V + BV,)- 
30 " A 30 


For body-source network. 


aw . A . 3 ft 1 3Q 

^ ^ ■ 30"^ 2 30 


(F.15) 


(F.16) 


3[-2(iL+ Bv'^)-((n®vu)®n)] 


“ -2B ■^•((ftOVli) ®n) -2 (V„+ 8 V.y 

d A " 3A 

—>■ 

3 V 

° + 2(C+ BT^)-(n®-^; 

3{-2( C-t BTfl)-((n®^lj)®lf)t , 0 -e.X.v 


(F.17) 

9((ft (E)^) (X) n) 


- A ^ 

-2B— •((n®VM)® n) - 2(Voo + B vI) 

30 A 
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II llll 


(5) 


3f 


ave* 




where 


1 ^ 

— n. 

Vie ^ 

1 ^ 
Me 


f. > 


f e* ^e”Ag ® ig 


-*■ 3?.\(F.19) 

^ (F.20) 

for f_ 


for fv 


(F.21) 


(6) and E is a summation overall fed sheet network panels. 

3v^.(Efi + fg) 


3f 

30 


0= 0,A,v 


30 


av. 


(F.22) 


30 


3f . 3f 

i 1' 


where V , f^ and fg are defined as above. 
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Appendix G 
Quasi-Newton Scheme 

Denote all geometry degrees of freedom by the vector 0 and all the 
singularity parameters by A . The boundary conditions determining the 
singularity strength are given by equation set F.l. Denote the set by 


F(A,0) - 0 


(G.l) 


Equations F.2 and F.3 determine the geometric degrees of freedom. Denote this 
set by 


G(A,0) - 0 


(G.2) 


Small perturbation of equations G.l and G.2 from the initial "starting 
solution" result in a set of linear equations governing the perturbation 
variables A , 0 . 



The perturbation quantities in equation G.3 are denoted symbolically as AX, 
the coefficient matrix (Jacobian) as J, and the right-hand side as -f . 
Equation (G.3) becomes 


JAX = -f 


(G.4) 


The set of equation G.l and G.2 is solved iteratively with a Quasi-Newton 
scheme. Represent the ith iteration by superscript i. The scheme proceeds to 
find the corrections AX(i) from the equation 


j(i) 


(G.5) 


and forms the new approximate solution (termed the next iterate) 


j^(i + 1) « x(i) + fiCi) 


(G.6) 
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where J(i) = J(X(i), f(i) = f(X(1)), and 6(1) is a scaling 

parameter to limit the step size of the correction vector. The Jacobian at 

x(l'^i) is obtained by using the following update formula (reference 38) 


j(i+l) jU)+dU) 


(G.7) 


where 


fi) T 

Ax(’)-rAx(^) ‘ ’ 

In this way, there is no need to reevaluate the partial derivatives comprising 
the elements of the Jacobian at each iteration. The superscript T denote the 
transpose of a vector. The aerodynamic influence coefficients which are 
changed by the geometry update are recalculated every iteration. 


The scaling parameter 6^^) is introduced to alleviate the problem of 
overshoot in the classical Newton scheme. For each iteration cycle, the 
following criteria are used to determine 6(^': 


0 < < 1 


(8.9) 


and 

6(i) 1^0(1) I < Y 


(G.IO) 


where Y is a predetermined quantity chosen to limit the maximum correction 
for the panel orientation angles (10 degrees in the code). In addition, a 
halving process of the scaling parameter 6(^) is applied to ensure the 
inequal ity 

II <11 (G.ii) 

where II II is the Euclidean norm representing the length of a vector. 


This halving process is performed repeatedly until either the criterion in 
equation 6.11 is satisfied or three cycles of the step size reduction (i.e., 
three halvings) are completed. The quality of the solution is monitored by 
examination of the sum of squares of residuals defined by 


114 



- IIFJ^ +11G||^ 


(G.12) 


To initiate the solution process, an initial geometry is required. The size 
of the fed sheet and the initial free sheet geometry are taken from Smith's 
conical flow solutions or, as experience allows, by assuming an initial 
geometry. 
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